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Abstract 

Anatomic connections between brain areas affect information flow between neuronal circuits 
and the synchronization of neuronal activity. However, such structural connectivity does not 
coincide with effective connectivity, related to the more elusive question "Which areas cause 
the present activity of which others?". Effective connectivity is directed and depends flexi- 
bly on contexts and tasks. Here we show that a dynamic effective connectivity can emerge 
from transitions in the collective organization of coherent neural activity. Integrating simula- 
tion and semi-analytic approaches, we study mesoscale network motifs of interacting cortical 
areas, modeled as large random networks of spiking neurons or as simple rate units. Through a 
causal analysis of time-series of model neural activity, we show that different dynamical states 
generated by a same structural connectivity motif correspond to distinct effective connectivity 
motifs. Such effective motifs can display a dominant directionality, due to spontaneous symme- 
try breaking and effective entrainment between local brain rhythms, although all connections in 
the considered structural motifs are reciprocal. We show then that transitions between effective 
connectivity configurations (like, for instance, reversal in the direction of inter-areal interac- 
tions) can be triggered reliably by brief perturbation inputs, properly timed with respect to 
an ongoing local oscillation, without the need for plastic synaptic changes. Finally, we analyze 
how the information encoded in spiking patterns of a local neuronal population is propagated 
across a fixed structural connectivity motif, demonstrating that changes in the active effective 
connectivity regulate both the efficiency and the directionality of information transfer. Previous 
studies stressed the role played by coherent oscillations in establishing efficient communication 
between distant areas. Going beyond these early proposals, we advance here that dynamic in- 
teractions between brain rhythms provide as well the basis for the self-organized control of this 
"communication-through-coherence" , making thus possible a fast "on-demand" reconfiguration 
of global information routing modalities. 

Author Summary 

The circuits of the brain must perform a daunting amount of functions. But how can "brain 
states" be flexibly controlled, given that anatomic inter-areal connections can be considered as 
fixed, on timescales relevant for behavior? We hypothesize that, thanks to the nonlinear in- 
teraction between brain rhythms, even a simple circuit involving few brain areas can originate 
a multitude of effective circuits, associated with alternative functions selectable "on demand". 
A distinction is usually made between structural connectivity, which describes actual synap- 
tic connections, and effective connectivity, quantifying, beyond correlation, directed inter-areal 
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causal influences. In our study, we measure effective connectivity based on time-series of neural 
activity generated by model inter-areal circuits. We find that "causality follows dynamics". 
We show indeed that different effective networks correspond to different dynamical states as- 
sociated to a same structural network (in particular, different phase- locking patterns between 
local neuronal oscillations). We then find that "information follows causality" (and thus, again, 
dynamics). We demonstrate that different effective networks give rise to alternative modalities 
of information routing between brain areas wired together in a fixed structural network. Thus, 
we propose that the flow of "digital-like" information encoded in spiking patterns is controlled 
by the self-organization of interacting "analog" rate oscillations. 



Introduction 

In Arcimboldo's (1527-1593) paintings, whimsical portraits emerge out of arrangements of flow- 
ers and vegetables. Only directing attention to details, the illusion of seeing a face is suppressed 
(Fig. [T]^-B). Our brain is indeed hardwired to detect facial features and a complex network of 
brain areas is devoted to face perception jl]. The capacity to detect faces in an Arcimboldo 
canvas may be lost when lesions impair the connectivity between these areas [2]. It is not con- 
ceivable, however, that, in a healthy subject, shifts between alternate perceptions are obtained 
by actual "plugging and unplugging" of synapses, as in a manual telephone switchboard. 

Brain functions — from vision [s] or motor preparation [4] up to memory [s], attention |6-8] 
or awareness jol — as well as their complex coordination jlOj require the control of inter-areal 



interactions on time-scales faster than synaptic changes 11,12 . In particular, strength and 



direction of causal influences between areas — described by the so-called effective connectivity 



(or, in a more restrictive sense, causal connectivity) [13-16 — must be reconfigurable even when 
the underlying structural (i.e. anatomic) connectivity is fixed. The ability to quickly reshape 
effective connectivity is a chief requirement for performance in a changing environment. Yet 
it is an open problem to understand which circuit mechanisms allow for achieving this ability. 
How can manifold effective connectivities — corresponding to different patterns of inter-areal 
interactions, or brain states |17| — result from a fixed structural connectivity? And how can 
effective connectivity be controlled without resorting to structural plasticity, leading to a flexible 
"on demand" selection of function? 

Several experimental and theoretical studies have suggested that multi- stability of neural 
circuits might underlie the switching between different perceptions or behaviors [T8] - |22| . In this 
view, transitions between many possible attractors of the neural dynamics would occur under 
the combined influence of structured "brain noise" (23] and of the bias exerted by sensory or 



cognitive driving 24-26 . Recent reports have more specifically highlighted how dynamic multi- 



stability can give rise to transitions between different oscillatory states of brain dynamics 271281. 



This is particularly relevant in this context, because long-range oscillatory coherence |12||29| — in 



particular in the gamma band of frequency (30-100 Hz) (29 -32 — is believed to play a central 
role in inter-areal communication. 

Ongoing local oscillatory activity modulates rhythmically neuronal excitability [33j. As 
a consequence, according to the influential communication-through-coherence hypothesis (311, 



neuronal groups oscillating in a suitable phase coherence relation — such to align their respective 
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"communication windows" — are likely to interact more efficiently than neuronal groups which 
are not synchronized. However, despite accumulating experimental evidence of communication- 
through-coherence mechanisms 34 - 38 and of their involvement in selective attention and top- 
down modulation |30^,39, 40j , a complete understanding of how inter-areal phase coherence can 
be flexibly regulated at the circuit level is still missing. In this study we go beyond earlier 
contributions, by showing that the self-organization properties of interacting brain rhythms lead 
spontaneously to the emergence of mechanisms for the robust and reliable control of inter-areal 
phase-relations and information routing. 

Through large-scale simulations of networks of spiking neurons and rigorous analysis of 
mean-field rate models, we model the oscillatory dynamics of generic brain circuits involving 
a small number of interacting areas {structural connectivity motifs at the mesoscopic scale). 



Following 41 , we extract then the effective connectivity associated to this simulated neural 
activity. In the framework of this study, we use a data driven rather than a model driven 
approach to effective connectivity (TgI (see also Discussion section), and we quantify causal 



influences in an operational sense, based on a statistical analysis of multivariate time-series of 



synthetic "LFP" signals. Our causality measure of choice is Transfer Entropy (TE) [42,43 . TE 



is based on information theory [44] (and therefore more general than causality measures based 



on regression 45 461), is "model- agnostic" and in principle capable of capturing arbitrary linear 



and nonlinear inter-areal interactions. 

Through our analyses, we first confirm the intuition that "causality follows dynamics". In- 
deed we show that our causal analysis based on TE is able to capture the complex multi-stable 
dynamics of the simulated neural activity. As a result, different effective connectivity motifs 
stem out of different dynamical states of the underlying structural connectivity motif (more 
specifically, different phase-locking patterns of coherent gamma oscillations). Transitions be- 
tween these effective connectivity motifs correspond to switchings between alternative dynamic 
attr actors. 

We show then that transitions can be reliably induced through brief transient perturba- 
tions properly timed with respect to the ongoing rhythms, due to the non-linear phase-response 
properties 48 of oscillating neuronal populations. Based on dynamics, this neurally-plausible 
mechanism for brain-state switching is metabolically more efficient than coordinated plastic 
changes of a large number of synapses, and is faster than neuromodulation |49| . 

Finally, we find that "information follows causality" (and, thus, again, dynamics). As a 
matter of fact, effective connectivity is measured in terms of time-series of "LFP-like" signals 
reflecting collective activity of population of neurons, while the information encoded in neuronal 
representations is carried by spiking activity. Therefore an effective connectivity analysis — even 
when based on TE — does not provide an actual description of information transmission in the 
sense of neural information processing and complementary analyses are required to investigate 
this aspect. Based on a general information theoretical perspective, which does not require 
specifying details of the used encoding [44| , we consider information encoded in spiking pat- 
terns 50-54 , rather than in modulations of the population flring rate. As a matter of fact, the 



spiking of individual neurons can be very irregular even when the collective rate oscillations are 



regular [55-58 . Therefore, even local rhythms in which the firing rate is modulated in a very 
stereotyped way, might correspond to irregular (highly entropic) sequences of codewords encod- 
ing information in a digital- like fashion (e.g. by the firing — "1" — or missed firing — "0" — 



4 



of specific spikes at a given cycle ^9]). In such a framework, oscillations would not directly 
represent information, but would rather act as a carrier of "data-packets" associated to spike 
patterns of synchronously active cell assemblies. By quantifying through a Mutual Information 
(MI) analysis the maximum amount of information encoded potentially in the spiking activity 
of a local area and by evaluating how much of this information is actually transferred to dis- 
tant interconnected areas, we demonstrate that different effective connectivity configurations 
correspond to different modalities of information routing. Therefore, the pathways along which 
information propagates can be reconfigured within the time of a few reference oscillation cycles, 
by switching to a different effective connectivity motif. 

Our results provide thus novel theoretical support to the hypothesis that dynamic effective 
connectivity stems from the self-organization of brain rhythmic activity. Going beyond pre- 
vious proposals, which stressed the importance of oscillations for feature binding ^60j or for 
efficient inter-areal "communication-through-coherence", we advance that the complex dynam- 
ics of interacting brain rhythms allow to implement reconfigurable routing of information in a 
self-organized manner and in a way reminiscent of a clocked device (in which digital-like spike 
pattern codewords are exchanged at each cycle of an analog rate oscillation). 



Results 

Models of interacting areas 

In order to model the neuronal activity of interacting areas, we use two different approaches, 
previously introduced in [61j. First, each area is modeled as a large network of thousands of 
excitatory and inhibitory spiking neurons, driven by uncorrelated noise representing background 
cortical input (network model). Recurrent synaptic connections are random and sparse. In these 
networks, local interactions are excitatory and inhibitory. A scheme of the network model for 
a local area is depicted in Fig. [2]A. (left). In agreement with experimental evidence that the 
recruitment of local interneuronal networks is necessary for obtaining coherent gamma cortical 



activity in vitro and in vivo 62 63 , the model develops synchronous oscillations (~ 50 Hz) 



when inhibition is strong, i.e. for a sufficiently large probability pi of inhibitory connection 



[55-58 64|. These fast oscillations are clearly visible in the average membrane potential (denoted 
in the following as "LFP"), an example trace of which is represented in Fig. [2]A. (bottom right). 
Despite the regularity of these collective rhythms, the ongoing neural activity is only sparsely 



synchronized. The spiking of individual neurons is indeed very irregular [55,57 and neurons 
do not fire an action potential at every oscillation cycle, as visible from the example spike 
trains represented in Fig. [2]A (top right). Structural network motifs involving N > 2 areas 
are constructed by allowing excitatory neurons to establish in addition long-range connections 
toward excitatory or inhibitory neurons in a distant target area (see a schematic representation 
of an A'^ = 2 structural connectivity motif in Fig. [2]C). The strength of inter-areal coupling is 
regulated by varying the probability pE of establishing an excitatory connection. 

In a second analytically more tractable approach, each area is described by a mean-field firing 
rate variable (rate model). The firing rate of a local population of neurons obeys the non- linear 
dynamical equation Q (see Methods). All incorporated interactions are delayed, accounting for 
axonal propagation and synaptic integration. Local interactions are dominantly inhibitory (with 
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coupling strength Kj < and delay D). Driving is provided by a constant external current. 
A cartoon of the rate model for a local area is depicted in Fig. (left). As in the network 
model, the firing rates undergo fast oscillations for strong inhibition {Kj < Ky ~ — '61]). 
An example firing rate trace is shown in Fig. [2jB (right). In order to build structural networks 
involving N >2 areas, different mean-field units are coupled together reciprocally by excitatory 
long range interactions with strength Ke > and delay D > D (see a schematic representation 
of an = 2 structural motif in Fig. [2p). Remarkably, the rate model and the network model 
display matching dynamical states |61| (see also later. Figures [sj |4] ancQ . More details on the 
network and the rate models are given in the Methods section and in the Supporting Text SI. 



Causality follows dynamics 

For simplicity, we study fully connected structural motifs involving a few areas (A^ = 2,3). 
Note however that our approach might be extended to other structural motifs |65| or even to 



larger-scale networks with more specific topologies 41 , 66 . 

In the simple structural motifs we consider, delays and strengths of local excitation and 
inhibition are homogeneous across different areas. Long-range inter-areal connections are as 
well isotropic, i.e. strengths and delays of inter-areal interactions are the same in all directions. 
Delay and strength of local and long-range connections can be changed parametrically, but only 
in a matching way for homologous connections, in such a way that the overall topology of the 
structural motif is left unchanged. As previously shown in [61] , different dynamical states — 
characterized by oscillations with different phase-locking relations and degrees of periodicity — 
can arise from these simple structural motif topologies. Changes in the strength of local inhibi- 
tion, of long-range excitation or of delays of local and long-range connections can lead to phase 
transitions between qualitatively distinct dynamical states. Interestingly, however, within broad 
ranges of parameters, multi-stabilities between dynamical states with different phase-locking 
patterns take place even for completely fixed interaction strengths and delays. 

We generate multivariate time-series of simulated "LFPs" in different dynamical states of 
our models and we calculate TEs for all the possible directed pairwise interactions. We show 
then that effective connectivities associated to different dynamical states are also different. The 
resulting effective connectivities can be depicted in diagrammatic form by drawing an arrow 
for each statistically significant causal interaction. The thickness of each arrow encodes the 
strength of the corresponding interaction. This graphical representation makes apparent, then, 
that effective connectivity motifs or, more briefly, effective motifs, with many different topologies 
emerge from structural motifs with a same fixed topology. Such effective motifs are organized 
into families. All the motifs within a same family correspond to dynamical states which are 
multi-stable for a given choice of parameters, while different families of motifs are obtained for 
different ranges of parameters leading to different ensembles of dynamical states. 

We analyze in detail, in Figures |3j |4] and [Sj three families of motifs arising for strong intra- 
areal inhibition and similarly small values of delays for local and long-range connections. We 
consider N = 2 (panels A and B) and N = 3 (panels C and D) structural motifs. Panels A 
and C show TEs for different directions of interaction, together with "LFPs" and example spike 
trains (from the network model), and rate traces (from matching dynamical states of the rate 
model). Panels B and D display motifs belonging to the corresponding effective motif families. 
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A first family of effective motifs occurs for weak inter-areal coupling. In this case, neuronal 
activity oscillates in a roughly periodic fashion (Fig. |3j^ and C, left sub-panel). When local 
inhibition is strong, the local oscillations generated within different areas lock in an out-of-phase 
fashion. It is therefore possible to identify a leader area whose oscillations lead in phase over 
the oscillation of laggard areas |6l| . In this family, causal interactions are statistically significant 
only for pairwise interactions proceeding from a phase-leading area to a phase-lagging area, as 
shown by the the box-plots of Fig. [3]A. and C (right sub-panel, see Discussion and Methods for a 
discussion of the threshold used for statistical significancy) . As commented more in detail in the 
Discussion section, the anisotropy of causal influences in leader-to-laggard and laggard-to-leader 
directions can be understood in terms of the communication-through-coherence theory. Indeed 
the longer latency from the oscillations of the laggard area to the oscillations of the leader area 
reduces the likelihood that rate fluctuations originated locally within a laggard area trigger 
correlated rate fluctuations within a leading area [35] (see also Discussion). Thus, out-of-phase 
lockings for weak inter-areal coupling give rise to a family of unidirectional driving effective 
motifs. In the case of = 2, causality is significant only in one of two possible directions 
(Fig. |3^), depending on which of the two areas assumes the role of leader. In the case of A^ = 3, 
it is possible to identify a "causal source" area and a "causal sink" area (see for an analogous 
terminology), such that no direct or indirect causal interactions in a backward sense from the 
sink area to the source area are statistically significant. Therefore, the unidirectional driving 
effective motif family for A^ = 3 contains six motifs (Fig. [3p), corresponding to all the possible 
combinations of source and sink areas. 

A second family of effective motifs occurs for intermediate inter-areal coupling. In this case, 
the periodicity of the "LFP" oscillations is disrupted by the emergence of large correlated fiuc- 
tuations in oscillation cycle amplitudes and durations. As a result, the phase-locking between 
"LFPs" becomes only approximate, even if it continues to be out-of-phase on average. The 
rhythm of the laggard area is now more irregular than the rhythm in the leader area. Laggard 
oscillation amplitudes and durations in fact fluctuate chaotically (Fig. |4]A and C, left sub-panel). 
Fluctuations in cycle length do occasionally shorten the laggard-to-leader latencies, enhancing 
non-linearly and transiently the influence of laggard areas on the leader activity. Correspond- 
ingly, TEs in leader-to-laggard directions continue to be larger, but TEs in laggard-to-leader 
directions are now also statistically signiflcant (Fig. [4]A and C, right sub-panel). The associated 
effective motifs are no more unidirectional, but continue to display a dominant direction or sense 
of rotation (Fig. and D). We refer to this family of effective motifs as to a family of leaky 
driving effective motifs (containing two motifs for N = 2 and six motifs for A^ = 3). 

Finally, a third family of effective motifs occurs for stronger inter-areal coupling. In this case 
the rhythms of all the areas become equally irregular, characterized by an analogous level of 
fluctuations in cycle and duration amplitudes. During brief transients, leader areas can still be 
identified, but these transients do not lead to a stable dynamic behavior and different areas in the 
structural motif continually exchange their leadership role (Fig. [sjA and C, left sub-panel). As a 
result of the instability of phase-leadership relations, only average TEs can be evaluated, yielding 
to equally large TE values for all pairwise directed interactions (Fig. [sj^ and C, right sub-panel). 
This results in a family containing a single mutual driving effective motif (Fig. [5^ and D). 

Further increases of the inter-areal coupling strength do not restore stable phase-locking 
relations and, consequently, do not lead to additional families of effective motifs. Note however 
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that the effective motif famihes explored in Figures [3j |4] and [5] are not the only one that can 
be generated by the considered fully symmetric structural motifs. Indeed other dynamical 
configurations exist. In particular, anti-phase locking (i.e. locking with phase-shifts of 180° for 
N = 2 and of 120° for A'" = 3) would become stable when assuming the same interaction delays 
and inter-areal coupling strengths of Figures [3j|4] and [5j but a weaker local inhibition. Assuming 
different interaction delays for local and long-range interactions, out-of-phase lockings continue 
to be very common, but in-phase and anti-phase locking can become stable even for strong 
local inhibition, within specific ranges of the ratio between local and long-range delays [61j . 
For = 3, in the case of general delays, more complex combinations can arise as well, like, 
for instance, states in which two areas oscillate in-phase, while a third is out-of-phase. In- 
phase locking between areas gives rise to identical TEs for all possible directed interactions, 
resulting in effective motifs without a dominant directionality. Anti-phase lockings for A'^ = 
2,3 give rise to relatively large inter-areal phase-shifts and, correspondingly, to weak inter- 
areal influences (at least in the case of weak inter-areal coupling), resulting in small TE levels 
which are not statistically significant (not shown). However, in the framework of this study, we 
focus exclusively on out-of-phase-locked dynamical states, because they are particularly relevant 
when trying to achieve a reconfigurable inter-areal routing of information (see later results and 
Discussion section). 

To conclude, we remark that absolute values of TE depend on specific parameter choices 
(notably, on time-lag and signal quantization, see Methods). However, the relative strengths 
of TE in different directions — and, therefore, the resulting topology of the associated effec- 
tive motifs — are rather robust against changes of these parameters. Robustness of causality 
estimation is analyzed more in detail in the Discussion section. 

Spontaneous symmetry breaking 

How can asymmetric causal influences emerge from a symmetric structural connectivity? A fun- 
damental dynamical mechanism involved in this phenomenon is known as spontaneous symmetry 
breaking. As shown in [gT] , for the case of the N = 2 structural motif, a phase transition occurs 
at a critical value of the strength of inter-areal inhibition. When local inhibition is stronger 
than this critical threshold, a phase-locked configuration in which the two areas oscillate in anti- 
phase loses its stability in favor of a pair of out-of-phase-locking configurations, which become 
concomitantly stable. The considered structural motif is symmetric, since it is left unchanged 
after a permutation of the two areas. However, while the anti-phase-locking configuration, sta- 
ble for weak local inhibition, share this permutation symmetry with the full system, this is no 
more true for the out-of-phase-locking configurations, stable for strong local inhibition. Note, 
nevertheless, that the configuration in which leader and laggard area are inverted is also a stable 
equilibrium, i.e. the complete set of stable equilibria continue to be symmetric, even if indi- 
vidual stable equilibria are not (leading thus to multi-stability). In general, one speaks about 
spontaneous symmetry breaking whenever a system with specific symmetry properties assumes 
dynamic configurations whose degree of symmetry is reduced with respect to the full symmetry 
of the system. The occurrence of symmetry breaking is the signature of a phase transition (of 
the second order j47j), which leads to the stabilization of states with reduced symmetry. 

The existence of a symmetry-breaking phase transition in the simple structural motifs we 
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analyze here (for simplicity, we consider the N = 2 case) can be proven analytically for the rate 
model, by deriving the function r(A(/)), which describes the temporal evolution of the phase-shift 



A(j) between two areas when they are weakly interacting 48 

d{A(t>) 



dt 



r(A0) (1) 



The function r(A(/)) for the rate model is shown in the left panel of Fig. [6^. Stable phase 
lockings are given by zeroes of r(A(/)) with negative slope crossing and are surrounded by basins 
of attraction (i.e. sets of configurations leading to a same equilibrium), whose boundaries are 
unstable in- and anti-phase lockings (Fig. [6]a). For the network model, a function r(A(/») with 
an analogous interpretation and a similar shape, shown in the right panel of Fig. [6^, can be 
extracted from simulations, based on a phase description of "LFP" time-series (see Methods and 
Supporting Figure SIA). The analogous distribution of the zero-crossings of r(A(/>) and r(A0) 
results in equivalent phase-locking behaviors for the rate and network models. Thus spontaneous 
symmetry breaking leads to multi-stability between alternative out-of-phase-lockings and to the 
emergence of unidirectional effective driving within a symmetric structural motif. 

Control of directed causality 

Because of multi-stability, transitions between effective motifs within a family can be triggered 
by transient perturbations, without need for structural changes. We theoretically determine 
conditions for such transitions to occur. The application of a pulse of current of small intensity h 
advances or delays the phase (f) of the ongoing local oscillation (see Supporting Figure SIB). This 
is true for rate oscillations of the mean-field rate model, but also for "LFP" oscillations reflecting 
rhythmic synchronization in the network model. In the latter case, the collective dynamics is 
perturbed by synchronously injecting pulse currents into all of the neurons within an area. The 
induced phase shift Scf) depends on the perturbation strength h but also on the phase (p at which 
the perturbation is applied. For the network model, this 5(j){(p] h) can be measured directly from 
numeric simulations of a perturbed dynamics (see Methods and right panel of Fig. [6p). For the 
rate model, the phase shift induced by an instantaneous phased perturbation can be described 
analytically in terms of the Phase Response Curve (PRC) Z{(f)) = §^ 48 (see Fig. 6|), left, and 



Supporting Text SI). After a pulse, the phase-shift between two areas is "kicked out" of the 
current equilibrium locking A0 and assumes a new transient value A</)* (solid paths in Fig.[6p), 
which, for weak perturbations and inter-areal coupling, reads: 

A^* = A(t) + 5cp{^;h) Acp + hZ{^)] (2) 

where the approximate equality between square brackets holds for the mean-field rate model. 
If A(/)* falls into the basin of attraction of a different phase-locking configuration than A0, the 
system will settle within few oscillation cycles into an effective connectivity motif with a different 
directionality (dashed green path in Fig. [ojC). Even relatively small perturbations can induce 
an actual transition, if applied in selected narrow phase intervals in which the induced S(p{(p; h) 
grows to large values. For most application phases, however, even relatively large perturbations 
fail to modify the effective driving direction (dashed red path in Fig. [6]C), because the induced 
perturbation 54){4>; h) is vanishingly small over large phase intervals (Fig. [6p). This is a robust 
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property, shared by the two (radicahy different) models we consider here and — we hypothesize — 
by any local circuit generating fast oscillations through a mechanism based on delayed mutual 
inhibition. As a consequence, for a given perturbation intensity, a successful switching to a 
different effective motif occurs only if the perturbation is applied within a specific phase interval, 
that can be determined analytically from the knowledge of T{A(j)) and of Z{<j)) for the rate 
model, or semi-analytically from the knowledge of r(A(/)) and S(f>{(j); h) (see Methods). Fig. [g^-F 
reports the fraction of simulated phased pulses that induced a change of effective directionality 
as a function of the phase of application of the perturbation. The phase intervals for successful 
switching predicted by the theory are highlighted in green. We performed simulations of the 
rate (Figs. [6^-F, left column) and of the network (Figs. [6^-F, middle column) models, for 
unidirectional (Figs. [6^) and leaky driving (Figs. [6^) effective motifs. Although our theory 
assumes small inter-areal coupling and is rigorous only for the rate model, the match between 
simulations and predictions is very good for both models and families of motifs. 

In Figs. [6^-F, we perturb the dynamics of the laggard area, but changes in directionality 
can also be achieved by perturbing the leader area (Supporting Figure S2). Note also that, in 
the network model, direction switchings can take place spontaneously, due to noisy background 
inputs. Such noise-induced transitions, however, occur typically on time-scales of the order of 
seconds, i.e. slow in terms of biologic function, because the phase range for successful switching 
induction is narrow. 



Effective entrainment 

A second non-linear dynamic mechanism underlying the sequence of effective motifs of Fig- 
ures [3] and [4] is effective entrainment. In this phenomenon, the complex dynamics of neural 
activity seems intriguingly to be dictated by effective rather than by structural connectivity. 

We consider as before a rate model of = 2 reciprocally connected areas (Fig. [2p). In order 
to properly characterize effective entrainment, we review the concept of bifurcation diagram 
[68| . As shown in |61|, when the inter-areal coupling Ke is increased, rate oscillations become 
gradually more complex (cfr. Fig. [7]A.), due to the onset of deterministic chaos (see also [69] 
for a similar mechanism in a more complex network model). For small Ke, oscillations are 
simply periodic (e.g. Ke = 4). Then, for intermediate Ke (e.g. Ke = 7), the peak amplitudes 
of the laggard area oscillation assume in alternation a small number of possible values (period 
doubling). Finally, for larger Ke (e.g. Ke = 8.5), the laggard peak amplitudes fluctuate in a 
random-like manner within a continuous range. This sequence of transitions can be visualized by 
plotting a dot for every observed value of the peak amplitudes of oscillation cycles, at different 
values of Ke- The accumulation of these dots traces an intricate branched structure, which 
constitutes the bifurcation diagram (Fig. [7^). 

Bifurcation diagrams for the leader and for the laggard area are plotted in Fig. [7^ (top 
panel, in orange and green color, respectively). We compare these bifurcation diagrams with 
the analogous diagrams constructed in the case of two unidirectionally coupled oscillating areas. 
Qualitatively similar bifurcation sequences are associated to the dynamics of the laggard area 
(bidirectional coupling) and of the driven area (unidirectional coupling. Fig. [7^, bottom panel, 
green color), for not too strong inter-areal couplings. In the case of unidirectional coupling, the 
peak amplitudes of the unperturbed driver area oscillations do not fluctuate at all. Therefore, 
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the corresponding bifurcation diagram is given by a constant line (Fig.jTp, bottom panel, orange 
color) . In the case of bidirectional coupling, the peak amplitudes of the leader area oscillations 
undergo fluctuations, but only with a tiny variance. Thus, the corresponding bifurcation diagram 
has still the appearance of a line, although now "thick" and curved (zooming would reveal 
bifurcating branches). Note that, for unidirectional coupling, the structural connectivity is 
explicitly asymmetric. The periodic forcing exerted by the driving area is then known to entrain 
the driven area into chaos f70|. Such direct entrainment is the dynamical cause of chaos. On the 
other hand, for bidirectional coupling, the structural connectivity is symmetric. However, due 
to spontaneous symmetry breaking, the resulting effective connectivity is asymmetric and the 
system behaves as if the leader area was a driver area, entraining the laggard area into chaos 
being only negligibly affected by its back-reaction. Such effective entrainment can be seen as 
an emergent dynamical cause of chaos. Thus, the dynamics of a symmetric structural motif 
with asymmetric effective connectivity and of a structural motif with a matching asymmetric 
topology are equivalent. 

For a sufficiently strong inter-areal coupling, symmetry in the dynamics of the bidirectional 
structural motif is suddenly restored [6lj, in correspondence with a transition to the mutual 
driving family of effective motifs (Figure [s]). As a result, in absence of symmetry breaking, 
effective driving cannot anymore take place. Thus, for a too strong inter-areal coupling, the 
emergent anisotropy of effective connectivity is lost, and, with it, the possibility of a dynamic 
control of effective connectivity (at least via the previously discussed strategies). 



Information follows causality 

Despite its name. Transfer Entropy is not directly related to a transfer of information in the 
sense of neuronal information processing. The TE from area X to area Y measures indeed 
just the degree to which the knowledge of the past "LFP" of X reduces the uncertainty about 
the future "LFP" of Y 43 [71 . As a matter of fact, however, the information stored in neural 



representations must be encoded in terms of spikes, independently from the neural code used. 
Therefore, it is important to understand to which extent an effective connectivity analysis based 
on "macroscopic" dynamics (i.e. TEs estimated from "LFPs") can pretend to describe actual 
"microscopic" information transmission (i.e. at the level of spiking correlations). 

In order to address this issue, we first introduce a framework in which to quantify the amount 
of information exchanged by interacting areas. In the case of our model, rate fluctuations could 
encode only a limited amount of information, since firing rate oscillations are rather stereotyped. 
On the other hand, a larger amount of information could be encoded based on spiking patterns, 
since the spiking activity of single neurons is very irregular and thus characterized by a large 



entropy 44,59 . As illustrated by Fig. [8jA, a code can be built, in which a "1" or a "0" symbol 
denote respectively firing or missed firing of a spike by a specific neuron at a given oscillation 
cycle. Based on such an encoding, the neural activity of a group of neurons is mapped to digital- 
like streams, "clocked" by the ongoing network rhythm, in which a different "word" is broadcast 
at each oscillation cycle. Note that we do not intend to claim that such a code is actually used in 
the brain. Nevertheless, we introduce it as a theoretical construct grounding a rigorous analysis 
of information transmission. 

We focus here on the fully symmetric structural motif of = 2 areas of Fig.[2p. We modify 
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the network model considered in the previous sections by embedding into it transmission lines 
(TLs), i.e. mono-directional fiber tracts dedicated to inter-areal communication (see Fig. |8^). 
In more detail, selected sub-populations of source excitatory neurons within each area establish 
synaptic contacts with matching target excitatory or inhibitory cells in the other area, in a one- 
to-one cell arrangement. Synapses in a TL are strengthened with respect to usual synapses, by 
multiplying their peak conductance by a multiplier Ktl (see Methods section). Such multiplier 
is selected to be large, but not too much, in order not to affect the phase-relations between 
the collective oscillations of the two areas. Indeed, selecting a too large Ktl would lead to 
an in-phase-locking configuration in which collective dynamics is enslaved to the synchronous 
activity of source and target populations. As analyzed in the Supporting Figure S3, a suitable 
tuning of Ktl ensures that source-to-target neuron communication is facilitated as much as 
possible, without disrupting the overall effective connectivity (associated to the unperturbed 
phase- locking pattern). Note that such TL synapses are here introduced as a heuristic device, 
allowing to maximize the potential capacity of inter-areal communication channels [44| . However, 
due to the occurrence of consistent spike-timing relations in out-of-phase locked populations, it 
might be that spike-timing-dependent plasticity ^72j lead to the gradual emergence of subsets 
of synapses with substantially enhanced weight [73| , which would play a role in inter-circuit 
communication very similar to TL synapses. 

The information transmission efficiency of each TL, for the case of different effective motifs, is 
then assessed by quantifying the Mutual Information (MI) [44[|59| between the "digitized" spike 
trains of pairs of source and target cells (see Methods). Since a source cell spikes on average every 
five or six oscillation cycles, the firing of a single neuron conveys H ~ 0.7 bits of information per 
oscillation cycle. MI normalized by the source entropy H indicates how much of this information 
reaches the target cell, a normalized MI equal to unity denoting lossless transmission. As shown 
by Fig. [8p-D, the communication efficiency of embedded TLs depends strongly on the active 
effective motif. In the case of unidirectional driving effective motifs (Fig. [sjC), communication 
is nearly optimal along the TL aligned with the effective connectivity. For the misaligned TL, 
however, no enhancement occurs with respect to control (i.e. pairs of connected cells not belong- 
ing to a TL). In the case of leaky driving effective motifs (Fig. [8p), communication efficiency is 
boosted for both TLs, but more for the TL aligned with the dominant effective direction. For 
both families of effective motifs, despite the strong anisotropy, the communication efficiencies 
of the two embedded TLs can be "swapped" within one or two oscillation cycles, by reversing 
the active effective connectivity through a suitable transient perturbation (see Fig. [6^-F). The 
considered N = 2 structural motif acts therefore as a "diode" through which information can 
propagate efficiently only in one (dynamically reconfigurable) direction determined by effective 
connectivity. 

Discussion 

Mechanisms for effective connectivity switching 

We have shown that a simple structural motif of interacting brain areas can give rise to multiple 
effective motifs with different directionality and strengths of effective connectivity, organized into 
different families. Such effective motifs correspond to distinct dynamical states of the underlying 
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structural motif. Beyond this, dynamic multi-stability makes the controlled switching between 
effective motifs within a same family possible without the need for any structural change. 

On the contrary, transitions between effective motifs belonging to different families (e.g. a 
transition from a unidirectional to a leaky driving motif) cannot take place without changes in 
the strength of the delay of inter-areal couplings, even if the overall topology of the underly- 
ing structural motif does not need to be modified. Each specific effective motif topology (i.e. 
motif family) is robust within broad ranges of synaptic conductances and latencies, however if 
parameters are set to be close to critical transition lines separating different dynamical regimes, 
transitions between different families might be triggered by moderate and unspecific parameter 
changes. This could be a potential role for neuromodulation, known to affect the net efficacy 
of excitatory transmission and whose effect on neural circuits can be modeled by coordinated 



changes in synaptic conductances 174,75 



Note that dynamical coordination of inter-areal interactions based on precisely-timed syn- 
chronous inputs would be compatible with experimental evidence of phase-coding [76[|8l] , in- 
dicating a functional role for the timing of spikes relative to ongoing brain rhythms (stimulus- 
locked |82^|83j as well as stimulus-induced or spontaneous |84j ) . Note also that the time of firing 



is potentially controllable with elevated precision 185 - 87 and has been found to depend on the 



phase of LFPs in local as well as in distant brain areas 37 . 

In general, control protocols different from the one proposed here might be implemented in 
the brain. For instance, phased pulses might be used as well to stabilize effective connectivity 
in the presence of stronger noise. Interestingly, the time periods framed by cycles of an ongoing 
oscillation can be sliced into distinct functional windows in which the application of the same 
perturbation produces different effects. 

Finally, in addition to "on demand" transitions, triggered by exogenous — sensory-driven — or 
endogenous — cognitive-driven — control signals, noise-driven switching between effective motifs 



might occur spontaneously, yielding complex patterns of activity during resting state 26,88,89 



Transfer Entropy as a measure of effective connectivity 

As revealed by our discussion of spontaneous symmetry breaking and effective entrainment, an 
effective connectivity analysis based on TE provides a description of complex inter-areal interac- 
tions compliant with a dynamical systems perspective. It provides, thus, an intuitive represen- 
tation of dynamical states that is in the same "space" as anatomical connectivity. Furthermore, 
as indicated by the analysis of Figure |8p-D, effective connectivity motifs obtained through TE 
point also at the specific modalities of information routing enabled by the associated dynamical 
states. In this sense, we can conclude that "causality follows dynamics". 

TE constitutes a model-free approach (although, non "parameter- free" , cfr. forthcoming 
section and Figure |9]) to the effective connectivity problem, suitable for exploratory data- 
driven analyses. In this sense it differs from regression-based methods like Granger Causality 



(GC) [45 ,46 or from Dynamic Causal Modeling (DCM) [90], which are model-driven 15 ,16 ,91 
Strategies like DCM, in particular, assume prior knowledge about the inputs to the system and 
works by comparing the likelihood of different a priori hypotheses about interaction structures. 
Such an approach has the undeniable advantage of providing a description of actual mechanisms 
underlying effective connectivity changes (the stimulus-dependence of effective couplings is in- 
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deed modeled phenomenologically) . However, it might be too restrictive (or arbitrary) when 
the required a-priori information is missing or highly uncertain. TE, on the contrary: does not 
require any hypothesis on the type of interaction; should be able to detect even purely non-linear 
interactions and should be robust against linear cross-talk between signals |92j. These features, 
together with the efficacy of TE for the causal analysis of synthetic time-series, advocate for a 
more widespread application of TE methods to real neural data [93| - |95| (at the moment limited 
by the need of very long time-series [92j ) . 

Note that we do not intend to claim superiority of TE in some general sense. As a matter of 



fact TE is equivalent to GC, as far as the statistics of the considered signals are gaussian 71 



Furthermore, non-linear generalizations of GC and DCM 96-99 might be able to capture at 
a certain extent the complex self-organized dynamics of the neural activity models analyzed in 
the present study. However, a systematic comparison of the performance of different methods 
in capturing causal connectivity of realistic non-linear models of neural dynamics goes beyond 
the focus of the present study and is deferred to future research. 

We finally would like to stress, to avoid any potential confusion, that the structural motifs 
analyzed in the present study are well distinct from causal graphical models of neural activity. 



in the statistical sense proper of DCMs 190 100 . They constitute indeed actual mechanistic 



models of interacting populations of spiking neurons, with a highly non-linear dynamics driven 
by background noise. Connections in these models are model synapses, i.e. mere structural cou- 
plings, not phenomenological effective couplings. Thus, effective connectivity is not constrained 
a priori, as in DCMs, but is an emergent property of network dynamics, consistent with the 
existence of effective motif topologies different from the underlying structural topology. 

Robustness of Transfer Entropy estimation 

The effective connectivity analyses presented in this study were conducted by evaluating TEs 
under specific parameter choices. However, absolute values of TE depend on parameters, like, 
notably, the resolution at which "LFP" signals are quantized and the time-lag at which we 
probe causal interactions. As discussed in detail in the Methods section, estimation of TE 
requires the sampling of joint distributions of "LFP" values in different areas at different times. 
Such distributions are sampled as histograms, based on discrete multi-dimensional binning. In 
practice, each "LFP" time-series is projected to a stream of symbols from a discrete alphabet. 



corresponding to different quantization levels of the continuous "LFP" signals 101 . The actual 
number B of used bins is a free parameter, although some guiding criteria for its selection do 
exist [43j. Concerning time-lag r, our TE analysis (conducted at the first Markov order [42j , 
following [4l||94]) describes predictability of "LFPs" at time t based on "LFPs" at time t — t. 
The used time-lag r is once again a free parameter. To deal with this arbitrariness in parameter 
choices, we explore systematically the dependence of TE estimations from the aforementioned 
parameters, by varying both B and r in a wide continuous range. Figure [9] summarizes the 
results of this analysis, for three different effective motifs. 

Considering the dependency on time-lag r, a periodic structure is clearly noticeable in the TE 
matrices reported in Figure [9j TE values tend to peak in precise bands of r, related to latencies 
between the oscillations of different areas. The analysis of the unidirectional driving motif 
(Figure [9]A) , associated to leader- laggard periodic configurations, is particularly transparent 
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(and has a high pedagogic value). Two characteristic time-lags can be defined: a "short" lag 
'TXY^ given by the time-shift from oscillation peaks of the leader area X to oscillation peaks of 
the laggard area Y; and a "long" lag, tyx = T — txy, given by the time-shift from laggard to 
leader oscillation peaks (here, T is an average oscillation period, common to both areas leader 
and laggard areas X and Y). TE in the direction from leader to laggard, TExy, peaks for the 
first time at a time-lag r = txy (and then at lags txy + nT, where n is a positive integer). 
TE in the direction from laggard to leader, TEyx-, peaks first at a time-lag r = Tyx (and 
then at lags Tyx + nT). If the "LFP" signals were deterministic and strictly periodic, the 
quantities TExYiTxY) and TEyxiryx) would be identical (and diverging for infinite precision 
[42| ) . However "LFP" signals are only periodic on average and have a stochastic component, 
due to the joint effect of random network connectivity and noisy background inputs. This 
stochastic component is responsible for small cycle-to-cycle fluctuations in the amplitude of 
"LFP" oscillation peaks. As discussed more in depth in a next subsection, the efficiency with 
which fluctuations in the output of a local area can induce (i.e., can "cause") fluctuations of 
the output of a distant interconnected area depends on the instantaneous local excitability 
of this target area, which is undergoing a rhythmic modulation due to the ongoing collective 



oscillation 31,33]. As a result, TE can reach different peak values in different directions (and, 
as a matter of fact, TExy{txy) > TEyxiryx))- 

Considering then the dependence on signal quantization, we observe that TE values tend to 
grow for increasing number of bins B, i.e. for a finer resolution in tracking "LFP" amplitude 
variations. This can be once again understood in terms of the temporal structure of "LFP" 
signals. As just mentioned, dynamic correlations between small "LFP" amplitude fluctuations 
carry information relevant for causality estimation. This information would be completely lost 
by using a too small number of bins for TE evaluation, given that the largest contribution to 
the dynamic range of "LFP" signals is provided by its fairly stereotyped oscillatory component. 
Conversely, using a too large number of bins would lead to under-sampling artifacts (therefore, 
we do not consider the use of more than B = 200 quantization bins). 

By evaluating a threshold for statistical significance independently for each direction and 
combination of B and r, we find that, for weak inter-areal coupling, TE never goes above this 
threshold in the laggard-to- leader direction (Figure |9j^). We are also unable to find any choice of 
B and r such that, for intermediate inter-areal coupling, TE in the laggard-to-leader direction 
becomes larger or equal than TE in leader-to- laggard direction (Figure |9^). Looking at matrices 
of the causal unbalancing ATE (see Methods, and Figure [9j third column), we see indeed that, 
for weak and intermediate coupling strengths, effective connectivity is robustly asymmetric in the 
parameter regions in which causal interactions are statistically significant. Effective connectivity 
is on the contrary balanced for strong inter-areal coupling (Figure [9]C). 

We can thus summarize the previous statements by saying that absolute values of TE depend 
on the choices of B and r, but that the topology of the resulting effective motif does not (at 
least in the wide range considered for this robustness analysis). 

Self-organized control of communication-through-coherence 

Traditionally, studies about communication-through-coherence or long-range binding between 
distant cell assemblies have emphasized the importance of in-phase locking (see, e.g. 35 , 102| ). 
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Although, as previously mentioned, in-phase locking (as well as anti-phase locking) can also arise 
in our models for different choices of coupling delays and inhibition strengths [6l] , we decided in 
the present study to focus on out-of-phase lockings. The case of spontaneous symmetry breaking 
is indeed particularly interesting, because it underlie the emergence of a dominant directionality 
in the causal influences between areas reciprocally coupled with comparable strengths. Further- 
more, spontaneous symmetry breaking is responsible for the multi-stability between effective 
connectivity configurations, thus opening the way to a self-organized control of inter-areal in- 
teractions 



In particular, our study confirms that the reorganization of oscillatory coherence might 
regulate the relative weight of bottom- up and top-down inter-areal influences [l7,30 or select 
different interaction modes within cortical networks involving areas of similar hierarchical level, 



as in the case of motor preparation or planning 4, 103 or language 104 



As a next step, we directly verified that "information follows causality", since changes in 
effective connectivity are paralleled by reconfiguration of inter-areal communication modalities. 
Following [32|35 , we explain the anisotropic modulations of communication efficiency (see Fig. [s]) 
in terms of a communication-through- coherence mechanism. In fact, because of the out-of-phase 
locking between rhythms, spikes emitted by neurons in a phase-leading area reach neurons in a 
phase-lagging area at a favorable phase in which they are highly excitable. Conversely, spikes 
emitted by neurons in a phase-lagging area reach neurons in a phase-leading area when they are 
strongly hyperpolarized by a preceding peak of synchronous inhibition. This same mechanism 
underlie also the anisotropy of "LFP"-based TE, since "LFP" fluctuations are the manifestation 
(at least in our model) of local population firing rate fluctuations. 

Therefore, by combining TE analyses of "LFP" -based effective connectivity with MI analyses 
of spike-based information transmission, we are able to establish a tight link between control 
of effective connectivity and control of communication-through-coherence, both of them being 
emergent manifestations of the self-organized dynamics of interacting brain rhythms. 

To conclude, we also note that similar mechanisms might be used beyond the mesoscale level 
addressed here. Multi-stabilities of structural motifs might be preserved when such motifs are 
interlaced as modules of a network at the whole- brain level 65 . Likewise, dynamic control of 



information routing between neuronal clusters 73 105 or even single cells might occur within 
more local microcircuits Il06l|107 . 



Communication-through-coherence beyond rate coding 

The previous discussions suggest that oscillations, rather than playing a direct role in the rep- 
resentation of information, would be instrumental to the reconfigurable routing of information 
encoded in spiking activity. Original formulations of the communication-through-coherence hy- 
pothesis suggested that oscillatory coherence facilitates the transmission of local fluctuations 
of firing rate to a distant site, thus assuming implicitly a rate-based encoding of information in 
neuronal activity. However, more complex coding mechanisms based on patterns of precisely 



timed spikes might be achievable by biologically-plausible neuronal circuits 85 , 86 



As a matter of fact, our study reveals that the inherent advantages of "labelled- line" codes 
[52[|108| (in which the information about which local neuron is firing is preserved) — i.e., notably, 
an augmented information capacity with respect to "summed-population" codes — might be 
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combined with the flexibiHty and the rehabiUty offered by the communication-through-coherence 
framework. Indeed, as shown by the analyses of Figure [8| suitable inter-areal phase relations 
make possible the transmission of information encoded in detailed spiking correlations, rather 
than just in population firing rate fluctuations. 

This is particularly interesting, since many cortical rhythms are only sparsely synchronized, 
with synchronous oscillations evident in LFP, Multi-Unit Activity or intracellular recordings but 
not in single unit spike trains 109 -111 . Such sparse firing might possibly reflect population- 



coding of behaviorally- relevant information transcending rate-based representations ^-54 . In- 
dependently from the complexity of these hypothetic representations, our study shows that self- 
organized communication-through-coherence would have the sufficient potential to dynamically 
route the rich information that these representations might convey. 



Perspectives 

It is very plausible that flexible inter-areal coordination is achieved in the brain through dy- 
namic self-organization 
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as in our models. However, qualitatively different mechanisms than 
symmetry breaking might contribute to the generation of dynamic effective connectivity in other 
regimes of activity. Despite sparse synchronization, the level of coherence in our model neuronal 
activity is larger than in many brain oscillations. However, our results might be generalized to 
activity regimes in which synchronization is weaker. Phase-relations have been shown to impact 



effective connectivity even in essentially asynchronous regimes 1121. It would be interesting to 



understand whether the dominant directionality of effective connectivity can be controlled when 
out-of-phase locking is only transient 



12,41 



Another open question is whether our theory can be extended to encompass the control of 

This is an important question since 



effective connectivity across multiple frequency bands 94 



top-down and bottom- up inter-areal communication might exploit different frequency channels, 
possibly due to different anatomic origins of ascending and descending cortico-cortical connec- 
tions 
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Finally, we are confident that our theory might inspire novel experiments, attempting to ma- 
nipulate the directionality of inter-areal influences via local stimulation applied conditionally to 
the phase of ongoing brain rhythms. Precisely timed perturbing inputs could indeed potentially 
be applied using techniques like electric 
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or optogenetic 
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in closed-loop implementations with millisecond precision 116 117] . 



microstimulation, especially 



Methods 

Network model 

Each area is represented by a random network of n^; 
Wang-Buzsaki-type conductance-based neurons 



4000 excitatory and nj = 4000 inhibitory 
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The Wang-Buzsaki model is described 
by a single compartment endowed with sodium and potassium currents. Note that results 
(not shown) of simulations performed with a more realistic ratio of n^; = 4000 excitatory and 
nj = 1000 inhibitory neurons per population would lead to qualitatively similar results with 
small parameter adjustments (using, for instance, parameters as in [69] ) . 
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The membrane potential is given by: 

dV 

C — = —II — iNa — Ik + hxt + Irec (3) 

where C is the capacitance of the neuron, II is a leakage current, lext is an external noisy driving 
current (due to background Poisson synaptic bombardment), and iNa and Ik are respectively a 
sodium and a potassium current, depending non linearly on voltage. The last input term I^ec is 
due to recurrent interactions with other neurons in the network. Excitatory synapses are of the 
AMPA-type and inhibitory synapses of the GABAA-type and are modeled as time-dependent 
conductances. A complete description of the model and a list of all its parameters are given in 
the Supporting Text SI. "LFP" A(t) = (V{t)) is defined as the average membrane potential 
over the Ne + Nj cells in each area. 

Short-range connections within a local area k from population Uk to population /3fc are 
established randomly with probability p^j^ , where a and /3 can be either one of the type E 
(excitatory) or I. The excitatory populations E^aie allowed also to establish connections toward 
populations Ei and Ii in remote areas {k ^ I). Such long-range connections are established with 
a probability (a = E,I). For simplicity, however, we assume that pfj = = pj and that 
P^EE ~ Pm ~ Pee ~ Pei — Pe- For each of the considered dynamical states, probabilities of 
connection are provided in the corresponding figure caption. 

Network model with embedded transmission lines (TLs) 

First, a structural motif of interconnected random networks of spiking neurons is generated, 
as in the previous section. Then, on top of the existing excitatory long-range connections, 
additional stronger long-range connections are introduced in order to form directed transmission 
lines. In each area a source sub-population, made out of 400 excitatory neurons, and a non- 
overlapping target sub-population, made out of 200 excitatory and 200 inhibitory neurons, are 
selected randomly. Excitatory cells in the source populations get connected to cells in the target 
sub-populations of the other area via strong synapses. These connections are established in a 
one-to-one arrangement (e.g. each source cell establishes a TL-synapse with a single target cell 
that does not receive on its turn any other TL-synapse). 

The peak conductance qtl of TL-synapses is Ktl times stronger than the normal excitatory 
peak conductance qe- For the simulations with TL (Fig. [sjof the main paper), we set Ktl = 
22 and 24.5 respectively for the unidirectional and for the leaky driving effective motifs. Such 
unrealistically strong peak conductances, whose purpose is to optimize information transfer by 
enhancing spiking correlations, can be justified by supposing that each source neuron establishes 
multiple weaker synaptic contacts with the same target neuron. The multiplier Ktl is selected 
to be as large as possible without altering the original out-of-phase locking relations between 
the two populations (Figure S3 A). Concretely, Ktl is tuned by raising it gradually until when 
a critical point is reached in which the populations lock in-phase (Figure S3C). Then, Ktl is 
set to be just below this critical point (Figure S3B). 
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Rate model 

Each area is represented by a single rate unit. The dynamical equations for the evolution of the 
average firing rate i?fe(t) in an area k are given by: 

Rkit) = -Rk{t) + [I + KiRk{t-D)+ (4) 
+ Y,KERiit-D)]+, k,l = l...N 

l^k 

Here, = x if X > 0, and zero otherwise. A constant current / represents a background 
input, Kj stands for the strength of intra-areal inhibition, Ke for the strength of inter-areal 
excitation and D and D are the delays of the local and long-range interactions, respectively. We 
consider in this study only fully symmetric structural motifs of mutually connected areas. 
For each of the considered dynamical states, the values of Kj, Ke, D and D are provided in the 
figure caption. 

Phase reduction and response 

Given an oscillatory time-series of neuronal activity, generated indifferently by a rate or by a 
network model, a phase = 360° • {t - tmax,i) / {tmax/+i - tmax/), for tmax,e <t< tmax,e+i, 
is linearly interpolated over each oscillation cycle. Here tmax/ denotes the start time of the 
^-th oscillation cycle. Note that this definition does not require that the oscillation is periodic: 
this empiric phase "elastically" adapts to fluctuations in the duration of oscillation cycles (see 
Supporting Figure SlA). 

The phase shift induced by a pulse perturbation I = h5{(p — (j)o) (see Supporting Figure SIB) 
is described by the Phase Response Curve (PRC) Z{(f>) = dcjy/dh (see Eq. ^ and [48]). For the 
rate model, the PRC can be evaluated analytically if certain general conditions on the relation 
between the oscillation period T and the local inhibition delay D are fulfilled |61j. Analytical 
expressions for the PRC of the rate model, as plotted in Fig. |6jD (left), are reported in the 
Supporting Text SI. 

In the network model, it is possible to evaluate the phase-shift induced by a perturbation, by 
directly simulating the effects of this perturbation on the oscillatory dynamics. A perturbation 
consists of a pulse current of strength h injected synchronously into all neurons of one area at a 
phase (f) of the ongoing local oscillation. The induced phase-shift (5(/>((/>; h) is estimated by com- 
paring the phases of the perturbed and of the unperturbed oscillations, when a new equilibrium 
is reached after the application of the perturbation. In detail, since the "LFP" time-series are 
not strictly periodic and the phase relation is fixed only on average, the average time-lag between 
the perturbed and the unperturbed "LFPs" is measured by computing their crosscorrelogram 
over 50 oscillation cycles, starting from the 10-th cycle after the perturbation. This average time 
lag (readable from the position of the crosscorrelogram peak) is then translated into a phase- 
shift, by dividing it by the average period (estimated through autocorrelation analysis of the 
perturbed and unperturbed time-series over the same observation window). Vanishingly small 
perturbations do not induce long-lasting phase-shifts. Therefore, moderately large perturbation 
strengths have to be used. In this case, the dependence of 6(j) on h is sensibly non-linear. As 
a consequence, we evaluate directly the resulting 5(j){(j)', h) for the used perturbation strength /i. 
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plotted in Fig. [6p (right). The quahtative shape of 5</>((/); h) however does not depend strongly 
on h. In particular, changes of h affect the amplitude of the maximum phase-shift but not the 
perturbation phase for which it occurs. The curve S(f){4>; h) is evaluated point- wise by applying 
perturbations at 100 different phases within a cycle. For each given phase, the perturbation 
is applied 100 times to 100 different cycles and the corresponding phase-shifts are averaged. 
Confidence intervals for 6<j3{(p; h) are determined phase-by-phase by finding the 2.5-th and the 
97.5-th percentile of the induced phase-shift distribution across these 100 trials. 

Phase locking 

For simplicity, we focus in the following on the case of = 2 areas, although our approach can 
be extended to larger motifs. The instantaneous phase-difference between two areas X and Y 
is given by A(/)(t) = mod \(j)x{t) — '?^y(i)) 360°]. For vanishing inter-areal coupling, the time 
evolution of Ai?i>(t) is described by Eq. ([T|. The term r(A(/)) is a functional of the phase response 
and of the limit cycle waveform of the uncoupled oscillating areas. For the rate model, r(A0) is 
determined from analytic expressions of Z{(f)) and of the rate oscillation limit cycle R{4>) (note 
that the dependence on t is replaced by a dependence on (j) after phase- reduction) for Ke = 0. 
It can be expressed as r(A(/)) = C(A0) — C(— Ac;/)), with: 

,■360° 

C{A(I))= Z{(j))R{(j) + A(j)- D)d(p (5) 

Jo° 

The resulting expression is reported in the Supporting Text SI and plotted in Fig. [6^ (left). 
Given Eq. ([T]), the phase shifts ±A(/)o between the two areas X and Y in stable phase- locked 
states correspond to top-down zero-crossings of the functional T{A(p) (i.e. zeroes with negative 
tangent slope, P' < 0). 

For the network model, the waveform of "LFP" oscillations can be determined through 
simulations. Since not all oscillation cycles are identical, the limit cycle waveform is averaged 
over 100 different cycles — as for the determination of 6(j){(j); h) — to yield an average limit cycle 
{A{(j))). Then, it is possible to evaluate a functional r(A(/<; h) = C{A(p; h) — (7(— A(/); /i), where: 

/>360° 

C{A(t>)= 5(0;/i)(A(0 + A(/.-L>))# (6) 

Jo° 

The functional T{A(j); h) is plotted in Fig. [g^S (right) for the used perturbation strength h. 
Although Eq. ([T]) does not exactly hold for the network model, the top-down zero-crossings 
of the functional r(A0; h) (whose position only weakly depends on h) continue to provide an 
approximation of the phase shifts ±A0o between the two areas X and Y in stable phase-locked 
states. In particular it is possible to predict whether the stable lockings will be in-phase, anti- 
phase or out-of-phase. 

Phase intervals for effective connectivity switching 

Phase intervals in which the application of a pulsed perturbation leads to a change of effective 
connectivity directionality are determined theoretically as shown below. For N = 2 and in a 
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given phase-locking state, the phase of the leader area can be written as 4>ieader = ^ + ^(po and 
the phase of the laggard area as (piaggard = ^- The application of a pulse perturbation of strength 
h to the laggard area shifts the phase of the ongoing local oscillation to 4>laggard ~ + '^(v^! ^)) 
where 5{ip;h) ~ hZ{ip) holds for the rate model in the case of small perturbations. If the 
achieved transient phase-shift between the two areas, Ac/)* ~ A(/)o — 5{<p;h), is falling into the 
basin of attraction of an alternative stable phase-locking (see Fig. [6p), then a switching toward a 
different effective motif takes place. Considering the dynamics of the instantaneous phase-shift, 
determined by the functionals T{A(j)) for the rate model and r(A(/>; h) for the network model 
(see Fig. [6^), switching will occur when: 

A0O < S{^p; h) < A(j)o + 180° (7) 

Here, we consider perturbations which induce a phase advancement, because the positive part 
of both the PRC in the rate model and the empiric 64>{(j); h) in the network model is larger than 
the negative part (see Fig. |6p). For a fixed perturbation intensity h, the condition ([T]) will be 
fulfilled only if when the phase ip of application of the perturbation falls within specific intervals, 
determined by the precise form of 6{(p; h). These intervals are highlighted in green in Fig. [6^ 
and F. Analogous considerations can be done in order to determine the intervals for successful 
switching when perturbing the leader area (see Supporting Figure S2). 

Transfer Entropy (TE) 

Let us consider first a structural motif with N = 2 areas. Let Ax(i) and Ay(t) be the 
"LFP" time-series of the two areas X and Y, and let quantize them into B discrete levels 
ii,...,iB (bins are equally sized). The continuous- valued "LFP" time-s eries are thus con- 
verted into strings of symbols Axit) and Ay(t) from a small alphabet 101 . Two transi- 
tion probability matrices, {PxY,Y{'^))ijk ~ ^i^vit + t) = £i\AY{t) = ij,Ax{t) = Ik] and 
(Py_y(r)).^. = P[Ay(t + r) = ^j|Ay(t) = where the lag r is an arbitrary temporal scale on 
which causal interactions are probed, are then sampled as normalized multi-dimensional his- 
tograms over very long symbolic sequences. These probabilities are sampled separately for each 
specific fixed phase-locking configuration. Epochs in which the system switches to a different 
phase-locking configuration, as well as transients following state switchings are dropped. The 
evaluation of Pxy,y{t) and Py^y(r) is thus based on disconnected symbolic subsequences, in- 



cluding overall O(IO^) oscillation cycles. Then, following 42 , the causal influence TExy(T) of 
area X on area Y is defined as the Transfer Entropy: 



TExy (r) = Pxy,y (r) log^ (8) 



where the sum runs over all the three indices i, j and k of the transition matrices. 

This quantity represents the Kullback-Leibler divergence [44| between the transition matri- 
ces Pxy,y{t) and Py_y(r), analogous to a distance between probability distributions. Therefore, 
TExy(T) will vanish if and only if PxY,YiT) and Py^y(T) coincide, i.e. if the transition proba- 
bilities between different "LFP" values of area Y do not depend on past "LFP" values of area 
X. Conversely, this quantity will be strictly positive if these two transition matrices differ, i.e. 
if the past "LFP" values of area X affect the evolution of the "LFP" in area Y. 
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We also measure the causal unbalancing 93 : 



ATE = ~ (9) 

which is normalized in the range — 1 < ATE < 1. A value close to zero denotes symmetric 
causal influences in the two directions, while large absolute values of ATE signal the emergence 
of asymmetric effective connectivity motifs. 

Partialized Transfer Entropy (pTE) 

Considering now a structural motif with = 3 areas, equation ([s]) has to be modified in order 
to distinguish causal interactions which are direct (e.g. X toward Y) from interactions which 
are indirect (e.g. X toward Y, but through Z). A solution allowing to assess only direct 



causal influences is partialization [42 ,71 . Indirect interactions from area X to area Y involving 
a third intermediate area Z are filtered out by conditioning the transition matrices for the 
"LFP" activity of Y with resepect to the activity of the Z. Two conditional transition matrices, 
{PxY,Y\z{r))^^^,, = P[Ay(t + r) = 4|Ay(t) = £„Ax{t) = £k,Az{t) = li] and (Py,y|z(T))^^, = 
P[Ay(t + r) = ^j|Ay(t) = ij, Az{t) = £i], are then constructed and used to evaluate: 

TExY\z{r) = J2 PxY,Y\z{r) logs ^^^^^^ (10) 

where the sum runs over all the four indices i, j, k and I. The effective connectivity in the panels 
C of Figures [sj |4] and [5] is computed using pTE according to equation (10). 



Statistic validation of effective connectivity 

Absolute values of TExy depend strongly on the time-lag r > and on the number of discrete 
levels B. Nevertheless, we find that relative strengths of causal influences are qualitatively un- 
changed over broad ranges of parameters, as displayed in the Supporting Figure SI. Furthermore 



the "plug-in" estimates of TE given by equations \m and ( 10 ) suffer from finite-sampling biases. 



and a rigorous debiasing procedure is not yet known 43 . Therefore, for each value of r and B it 



is necessary to assess the significancy of the inferred causal interactions through comparison with 



suitably randomly resampled data 119 . To estimate the confidence intervals for the estimated 



TEs and the baseline for significancy we adopt a geometric bootstrap method 120 , guaranteed 



to generate resampled time-series with similar auto- and cross-correlation properties up to a 
certain lag. This is important, since "LFP" time-series have a strong oscillatory component, 
whose correlation structure has to be maintained under resampling. Each resampled time-series 
A^(t) consists of a concatenation of blocks sampled from the original time-series Ax{t). Each 
A5|(t) has the same length as the original Ax{t). Every upward crossing, i.e. every time at 
which Ax{t) crosses from below its time-averaged value Ax(t), is a potential start-time for a 
block. The first element of each block is obtained by selecting randomly one of these potential 
start-times. Then, the block consists of the L oscillation cycles following the chosen start-time, 
where the random integer L follows a geometric distribution P{L) oc (1 — q)^~^, with < g < 1 
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and an average block length of {L) = 1/q (we have taken (L) = 20 oscillation cycles, longer than 
the mean correlation time for all the simulated "LFPs"). Randomly selected blocks are then 
concatenated up to the desired length. 

When considering a structural motif involving more areas, the "LFP" time-series of each area 
can be resampled jointly or independently. When resampling jointly, matching starting points 
and block-lengths are selected for each block of the resampled time-series of each area, leading 
to resampled multivariate time-series in which the structure of causal influences should not be 
altered. The distribution of TExy over jointly resampled "LFP" time-series describes then for 
each directed pair of areas X and Y the strength of the corresponding effective connectivity 
link, as well as the fluctuations of this strength. Conversely, when resampling independently 
the time-series, start-points and block-lengths of the resampled blocks are chosen independently 
for each area. This second procedure leads by construction to causally independent time-series. 
Any residual TE between directed pairs of independently resampled "LFPs" indicates therefore 
systematic biases, rather than actual causal influences. For each directed pairs of areas X 
and y, significance of the corresponding causal interaction can be assessed by comparing the 
bootstrapped distributions of TExy(T) and of TExy(T). This comparison is performed in 
Figures [3j|4] and [5] and in Supporting Figure S3D-E. Here, boxes indicate the median strength of 
TExy(T) for different directions and the corresponding confidence intervals, comprised between 
a lower extreme Qi — 1.5(Q3 — Qi) and and upper extreme + l.b{Q^ — Qi), where Qi, Q2 and 
Q3 are respectively the first, the second and the third quartiles of the distribution of TExy('7") 
over jointly resampled time-series. Median values of TE^y (t) and the corresponding confidence 
intervals, evaluated as before, are represented by horizontal dashed lines and a surrounding 
shaded band. When the distributions of TExy(T) and TEyx(T) are not significantly different, 
a single baseline band is plotted. In this study, strengths and base-line for significancy of effective 
connectivity for each direction are validated based on, respectively, 500 jointly resampled and 
500 independently resampled replicas. 

Note that geometric bootstrap can be applied to arbitrary signals, and does not depend on 
their strict periodicity. However it is precisely the strong periodic component of our signals that 
makes necessary the use of geometric bootstrap techniques. Indeed, conventional bootstrap, 
strongly disrupting signal periodicity, would lead to artificially low thresholds for statistical 
significance of TE (not shown) . 

Entropy and Mutual Information (MI) 

We evaluate information transmission between pairs of mono-synaptically connected cells in 
difi^erent areas, linked by a TL-synapse (TL pairs) or by a normally weak long-range synapse 
(control pairs). Inspired by |59| , spike trains are digitized into binary streams Sj(/c), where Sj(/c) 
= 1 or respectively when neuron i fires or does not fire during the k-th local oscillation cycle 
(cycle counting is performed independently for each area and includes all the oscillation cycles 
following a common reference initial time). Note that neurons fire very sparsely and, due to the 
elevated degree of synchrony in our model, only in narrow temporal intervals centered around the 
peaks of the ongoing "LFP" oscillations. In particular, they fire at maximum once per oscillation 
cycle. Thus, this oscillatory spiking activity is naturally quantized in time and binning [59j is 
not required. For each considered directed pair of cells {i source cell, j target cell), based on very 
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long duration spike trains, we sample normalized histograms for three probability distributions: 
Pi = P{si{k)), Pj = P{sj{k)) and Pij = P{si{k),Sj{k')). When sampling the joint probability- 
distribution Pij we have to distinguish two cases: (i) If the presynaptic cell i belongs to a leader 
area, i.e. the oscillation of the source area leads in phase over the oscillation of the target area of 
the considered synapse, then k' = k; (ii) Conversely, if the presynaptic cell i belongs to a laggard 
area, i.e. the oscillation of the target area leads in phase over the oscillation of the source area 
of the considered synapse, then k' = k + 1. This means that we seek for spiking correlations 
only in pairs of spiking (or missed spiking) events in which the "effect" follows temporally its 
potential "cause", since physical information transmission cannot occur backward in time. As 
for the estimation of TE (see previous section), the probabilities Pj, Pj and Pij are sampled 
separately for each specific phase- locking configuration of the ongoing "LFPs" . Epochs in which 
the system switches to a different phase-locking configuration, as well as transients following 
state switchings are dropped. The evaluation of these probabilities is thus based on disconnected 
spike train chunks, including overall O(IO^) oscillation cycles. Based on these probabilities, the 
Shannon entropy H of the spike train of the presynaptic neuron i (measuring the information 
content in its activity) is evaluated as: 

Hi = -^paog2Pi (11) 

and MI between pre- and postsynaptic cells as: 

MI is then normalized by the entropy of the pre-synaptic cell, in order to measure the relative 
efficiency of information transmission along each TL or control synapse. 

Statistics are taken over 400 pairs of cells per synapse set, i.e. one set of strong synapses per 
embedded TL, plus one set of (control) weak synapses. The box- plots in Fig.jsjC-D report median 
efficiencies of information transmission efficiencies (for different active effective connectivities), 
as well as their confidence intervals, estimated non-parametrically from distribution quartiles, 
as discussed above for TE. Both MI and H are computed for (finite) spike trains of the largest 
available length L. Following f 59|121 , it is possible to correct these results for finite-size sampling 



bias (see Supporting Figure S4). MI and H are computed again, based on randomly selected 
shorter matching sections of the full length spike trains. The results of MI / H obtained for 
various shorter lengths L/q are then plotted against the so-called inverse data fraction q, where 
g = 1 correspond then to estimations based on full length spike trains. Quadratic extrapolation 
to (7 = provides a debiased estimation of MI / H. Note that, in order to allow a more direct 
comparison with the non-debiased TE analysis, the results plotted in Fig. |8p-D do not include 
any finite-size correction. As a matter of fact, as discussed in Supporting Figure S4, finite size 
bias induces a small quantitative overestimation of information transmission efficiency (from 
~ 3% to ~ 6%), that does not affect qualitatively any of the results presented here. 
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Figure Legends 



Figure 1. Flexibility of brain function requires dynamic effective connectivity. This 
is illustrated by the example of a Giuseppe Arcimboldo's painting ( Vertumnus; 1590, Skoklosters 
Slott, Sweden). A: the illusion of seeing a face is due to the default activation of a network of 
brain areas dedicated to face recognition. B: however, selective attention to individual com- 
ponents — e.g. to a pear or a flower — suppresses this illusion by modulating the interaction 
between these and other brain areas. Therefore, effective connectivity, i.e. the specific active 
pattern of inter-areal influences, needs to be rewired "on demand" in a fast and reliable way, 
without changes in the underlying structural connectivity between the involved areas. 
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Figure 2. Models of interacting areas. A: in the network model, each local area is modeled 
as a large network of randomly and sparsely interconnected excitatory and inhibitory spiking 
neurons (inhibitory cells and synapses are in blue, excitatory cells and synapses are in red, ue = 
nj = 4000). Individual neurons spike irregularly (see the spike trains of eight representative 
neurons, top right), but the activity of the network undergoes a collective fast oscillation, visible 
in the average membrane potential (see example "LFP" trace, bottom right). B: in the rate 
model, each local area is modeled by a single mean-field rate unit with delayed local inhibition 
(of strength Kj < 0). Its dynamics, describing the average area activity, also undergoes a 
fast oscillation (see example rate trace, right). C-D: the interaction between multiple local 
areas (A'^ = 2 in the case of the reported graphical illustrations, green and orange shading 
indicate separate areas) is modeled by the dynamics: of multiple local spiking networks, mutually 
interconnected by long-range excitatory synapses (see panel C); or of multiple rate units, coupled 
reciprocally by delayed excitation (of strength Ke > 0, see panel D). 

Figure 3. Effective motifs of the unidirectional driving family. For weak inter-arcal 
coupling strengths, out-of-phase lockings of local periodic oscillations give rise to a family of 
"unidirectional driving" effective motif. The figure shows dynamics and corresponding effective 
connectivities for fully symmetric structural motifs with N = 2 (panels A-B) or A^ = 3 (panels C- 
D) areas. A: the dynamics of A" = 2 interacting areas (green and orange colors) is illustrated by 
"LFPs" (left, top row) and representative spike trains (left, middle row, two cells per each area) 
from the network model (horizontal bar is 20 ms, vertical bar is 20 mV), as well as by matching 
rate traces (left, bottom row) from the rate model (arbitrary time units). The right sub-panel 
reports the associated effective connectivity measured by Transfer Entropy (TE), evaluated from 
"LFPs" time-series, for all possible directed interactions (indicated by colored arrows). Boxes 
indicate the interquartile range and whiskers the confidence interval for the estimated TEs. TEs 
above the grey horizontal band indicate statistically significant causal influences (see Methods). 
B: to the right of the corresponding box-plot, effective connectivity is also represented in a 
diagrammatic form. Arrow thicknesses encode the strength of corresponding causal interactions 
(if statistically significant). Below this effective motif, a second motif in the same unidirectional 
driving family is plotted (with a smaller size), corresponding to another motif version with 
equivalent overall topology but reversed directionality. The parameters used for AT = 2 are, 
for the network model: pj = 0.25, pE = 0.01; and for the rate model: Kj = —250, Ke = 5, 
D = D = 0.1. C: this panels reports similar quantities as panel A, but now for a structural motif 
with N = S areas (green, orange and light blue colors). Effective connectivity is now measured 
by partialized Transfer Entropy (pTE; see Methods), in order to account only for direct causal 
interactions. D: the six effective motifs of the unidirectional driving family for N = 3 are also 
reported. The parameters used for A^ = 3 are, for the network model: pi = 0.33, pE = 0.006; 
and for the rate model: Kj = -300, Ke = 5, D = D = 0.1. 

Figure 4. Effective motifs of the leaky driving family. The figure shows dynamics and 
corresponding effective connectivities for fully symmetric structural motifs with N = 2 (panels 
A-B) or AT = 3 (panels C-D) areas, for intermediate inter-areal coupling strength, leading to 
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asymmetrically irregular oscillations, phase-locked with an average out-of-phase relation. A: 
the dynamics of = 2 interacting areas (green and orange colors) is illustrated by "LFPs" 
(left, top row) and representative spike trains (left, middle row, two cells per each area) from 
the network model (horizontal bar is 20 ms, vertical bar is 20 mV), as well as by matching 
rate traces (left, bottom row) from the rate model (arbitrary time units). The right sub-panel 
reports the associated effective connectivity measured by Transfer Entropy (TE), evaluated from 
"LFPs" time-series, for all possible directed interactions (indicated by colored arrows). Boxes 
indicate the interquartile range and whiskers the confidence interval for the estimated TEs. TEs 
above the grey horizontal band indicate statistically significant causal influences (see Methods). 
B: to the right of the corresponding box-plot, effective connectivity is also represented in a 
diagrammatic form. Arrow thicknesses encode the strength of corresponding causal interactions 
(if statistically significant). Below this effective motif, a second motif in the same unidirectional 
driving family is plotted (with a smaller size), corresponding to another motif version with 
equivalent overall topology but reversed directionality. The parameters used for N = 2 are, 
for the network model: pi = 0.25, pE = 0.09; and for the rate model: Kj = —250, Ke = 25, 
D = D = 0.1. C: this panels reports similar quantities as panel A, but now for a structural 
motif with N = 3 areas (green, orange and light blue colors) . Effective connectivity is measured 
by partialized Transfer Entropy (pTE; see Methods), in order to account for direct but not for 
indirect causal interactions. D: the six effective motifs of the unidirectional driving family for 
N = 3 are also reported. The parameters used for AT = 3 are, for the network model: pj = 0.33, 
PE = 0.06; and for the rate model: Kj = -300, Ke = 11, D = D = 0.1. 

Figure 5. Effective motifs of the mutual driving family. The figure shows dynamics and 

corresponding effective connectivities for fully symmetric structural motifs with N = 2 (panels A- 
B) or A^ = 3 (panels C-D) areas, for large inter-areal coupling strength, leading to symmetrically 
irregular oscillations, without a stable phase relation. A: the dynamics of A/" = 2 interacting 
areas (green and orange colors) is illustrated by "LFPs" (left, top row) and representative spike 
trains (left, middle row, two cells per each area) from the network model (horizontal bar is 20 
ms, vertical bar is 20 mV), as well as by matching rate traces (left, bottom row) from the rate 
model (arbitrary time units). The right sub-panel reports the associated effective connectivity 
measured by Transfer Entropy (TE), evaluated from "LFPs" time-series, for all possible directed 
interactions (indicated by colored arrows). Boxes indicate the interquartile range and whiskers 
the confidence interval for the estimated TEs. TEs above the grey horizontal band indicate 
statistically significant causal influences (see Methods). B: to the right of the corresponding 
box-plot, effective connectivity is also represented in a diagrammatic form. Arrow thicknesses 
encode the strength of corresponding causal interactions (if statistically significant). A single 
motif is included in this family The parameters used for N = 2 are, for the network model: 
PI = 0.25, PE = 0.15; and for the rate model: Kj = -250, Ke = 27, D = D = 0.1. C: 
this panels reports similar quantities as panel A, but now for a structural motif with N = 3 
areas (green, orange and light blue colors). Effective connectivity is measured by partialized 
Transfer Entropy (pTE; see Methods), in order to account for direct but not for indirect causal 
interactions. D: the mutual driving effective motif for AT = 3 is also reported. The parameters 
used for A^ = 3 are, for the network model: pi = 0.33, pE = 0.1; and for the rate model: 
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Ki = -300, Ke = 15, D = D = 0.1. 

Figure 6. Dynamic control of effective connectivity. A: symmetric structural motifs can 
give rise to asymmetric dynamics in which one area leads in phase over the other (spontaneous 
symmetry breaking). Basins of attraction (in phase-shift space) of distinct phase- locking con- 
figurations are schematically shown here (for N = 2). Empty circles stand for unstable in- and 
anti-phase lockings and filled circles for stable out-of-phase lockings (corresponding to unidi- 
rectional driving effective motifs). B: phase-shift evolution function T{A(j)) for the rate model 
(left, analytical solution, Kj = —250) and for the network model (right, numerical evaluation. 
Pi = 0.25). Empty and filled circles denote the same stable and unstable phase- lockings as in 
panel A. C: cartoon of successful (dashed green arrow) and unsuccessful (dashed grey arrow) 
switchings induced by brief perturbations (lightning icon). An input pulse to the system destabi- 
lizes transiently the current phase-locking (solid red and green arrows). For most perturbations, 
the system does not leave the current basin of attraction and the previous effective motif is 
restored (dashed red arrow). However, suitable perturbations can lead the system to switch 
to a different effective motif (dashed green arrow). D: a pulse of strength h induces a phase 
advancement of the collective oscillations, depending on its application phase cp, as described 
by the Phase Response Curve Z{(j)) (left, rate model; analytical solution, Kj = —250) or by 
the induced shift 5(f){(j);h) (right, network model; numerical evaluation, pj = 0.25). E-F: fre- 
quency histogram of successful switching for pulses applied at different phases (the laggard area 
is perturbed; h = 0.21 for the rate model and h = 500 pA for the network model). Predicted 
intervals for successful switching are marked in green, for the unidirectional (panel E) and for 
the leaky effective driving (panel F) motifs (left, rate model; right, network model; parameters 
as in Figures [3]and|4j). Dia grams of the induced transitions are shown in the third column (see 
SI, Figure S2 for perturbations applied to the leader area). 

Figure 7. Effective entrainment. A: examples of rate oscillations for different values of 
the inter-areal coupling in the rate model {Kj = —250, Ke = 4,7 and 8.5, from bottom to 
top). Filled circles denote peaks of oscillation cycles, different color fillings denote different 
peak amplitudes. B: The oscillatory dynamics is qualitatively altered by increasing inter-areal 
coupling, as visualized by bifurcation diagrams, constructed by plotting different peak ampli- 
tudes at constant Ke, as different dots (the dots corresponding to the peak amplitudes in panel 
A, are highlighted also here by filled circles of matching colors). Varying Ke in a continuous 
range, these dots trace a complex branched structure, denoting emergence of novel dynamical 
states. The bifurcation diagrams for the case of two symmetrically connected areas (top) and two 
unidirectionally connected areas (bottom) are very similar. For a symmetric structural motif, 
spontaneous symmetry breaking leads to effective entrainment, mimicking the direct entrain- 
ment, which occurs for an asymmetric unidirectional structural motif. Leader and laggard areas 
in effective entrainment behave similarly to the driver and driven area in direct entrainment 
(orange and green bifurcation diagrams, respectively). Note that different structural motifs give 
rise to equivalent effective motifs (see side diagrams). Note: a different version of panel B was 
previously published in [6l] as Supplementary Figure F. 
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Figure 8. Effective connectivity affects information propagation. A: in the case of 
sparsely synchronized oscillations, individual neurons fire irregularly (see four example spike 
trains, middle row) even when the local area activity undergoes a very regular collective rhythm 
(evident in "LFP" traces, bottom row). Therefore, a large amount of information can be po- 
tentially encoded, at every (analog) oscillation cycle, in the form of (digital-like) codewords in 
which "1" or "0" entries denote respectively firing or missed firing of a specific neuron in the 
considered cycle (top row). B: the strength of specific subsets of long-range excitatory synapses 
is systematically enhanced in order to form unidirectional "transmission lines" (TLs) embedded 
into the N = 2 symmetric structural motif (see Methods). Cells and synapses belonging to TLs 
are highlighted by pale green ( "green-to-orange" area direction) and lilac ( "orange-to-green" area 
direction) colors. Communication efficiency along TLs is quantified by the Mutual Information 
(MI) between spike trains of pairs of source and target cells connected directly by a TL synapse, 
normalized by the entropy (H) of the source cell. C-D: boxplots (see Figures § [4] and § of 
MI /H for different groups of interconnected cells and for different active effective motifs. Pale 
green and lilac arrows below the boxplots indicate pairs of cells interconnected by the TL marked 
with the corresponding color. A dot indicates control pairs of cells interconnected by ordinary 
weak long-range synapses. Green and orange arrows indicate the dominant directionality of the 
active effective connectivity motif. C: unidirectional driving effective motif family. Communica- 
tion efficiency is enhanced only along the TL aligned to the directionality of the active effective 
connectivity, while it is undistinguishable from control along the other TL. D: leaky driving 
effective motif family. Communication efficiency is enhanced along both TLs, but more along 
the TL aligned to the dominant directionality of the active effective connectivity. 

Figure 9. Transfer entropy depends on time lag and quantization. A-C: The matrices in 
these panels illustrate the dependence of TE (network model, N = 2 fully symmetric structural 
motif, cfr. Figures [sj |4] and [s]) on the number B of discretization bins used to describe the 
time-series of neural activity and on the adopted time lag Tiag between the time-series (see 
Methods). The matrices in the first two columns (from the left) report TEs in the two possible 
interaction directions, TExy and TEyx; and the matrices in the third column visualize the 
causal unbalancing ATE (—1 < ATE < 1), which quantifies the asymmetry between causal 
influences in the two directions (see Methods). All of these quantities are evaluated for different 
combinations of B and Ti^g. The vertical axes of the matrices correspond to the range 2 < 
B < 200 bins and the horizontal axes to the range 1 ms < Tiag < 60 ms. This range of time 
lags corresponds approximately to three oscillation periods. Horizontal scale lines indicate the 
average oscillation period ((T) = 16.4, 18.9 and 19.1 ms, respectively for panels A, B and C). 
Values of TE and ATE are color-coded (see color bars at the bottom, note the two different 
color scales for TE and ATE) . Black dotted lines in the matrices enclose regions in which TExy 
or TEyx rise above the threshold for significancy of the corresponding causal interaction (see 
Methods). These significance contours are overlayed in the corresponding ATE matrix. A star 
denotes the combination of B and Tiag used for the analysis throughout the main article (ri^g = 5 
ms, B = 175). Different rows report TE matrices for different effective motifs. A: unidirectional 
driving effective motif. B: leaky driving effective motif. C: mutual driving effective motif. 
Diagrams of these effective motifs are drawn in the fourth column as a visual reference. All 
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other parameters are as for the analyses of Figures [3j |4] and [5} 

Supporting Information Legends 

Supporting Figure SI. Phase reduction and phase response. A: oscillating time-series 
(in the example, a "LFP" time-series from the network model) can be described in terms of phase, 
even if they are not periodic in strict sense, by interpolating linearly an instantaneous empiric 
phase variable (p to the oscillation cycles (generally of unequal lengths). B: the application of 
a pulse current 51 induces a shift 5(j) in the oscillation phase of the ongoing oscillation (in the 
example, a rate trace from the rate model). The amplitude of the induced shift depends on the 
phase (p of the ongoing oscillation at which the perturbation is applied. 

Supporting Figure S2. Dynamic control of effective connectivity (perturbation ap- 
pUed to the leader area). A-B: frequency histogram of successful switching for pulses applied 
at different phases {h = 0.2/ for the rate model and h = 500 pA for the network model). Pre- 
dicted intervals for successful switching are marked in green, for the unidirectional (panel E) 
and for the leaky effective driving (panel F) motifs (left, rate model; right, network model; 
parameters as in Figures [3]and|4]). Dia grams of the induced transitions are shown in the third 
column (see Figure [6] for perturbations applied to the laggard area) . 

Supporting Figure S3. Effective connectivity vi^ith transmission lines (TLs). We 

consider a fully symmetric structural motif of = 2 structurally connected areas with embed- 
ded unidirectional TLs. Synapses involved in TLs are enhanced by multiplying the ordinary 
excitatory peak conductance by a multiplier Ktl- Raster plots relative to the spiking activity 
of excitatory neurons of the two areas are shown in panels A~C (green and orange color denote 
spikes of excitatory neurons from different populations, the horizontal scale line corresponds to 
20 ms) for a weak inter-areal coupling (unidirectional driving effective motif, see Figure |3] for 
parameters). A: when Kj-l = (no TL embedded), the synchronous oscillations of the two 
populations lock in an out-of-phase fashion. B: for Ktl = 22 (just below a critical value), the 
raster plot of the spiking activity is virtually indistinguishable from the raster plot of panel A. 
C: for Ktl = 22.5 (just above a critical value), the oscillations of the two populations lock in an 
in-phase configuration. D-E: Effective connectivities associated to different dynamical states are 
measured by Transfer Entropy (TE), evaluated from "LFPs" time-series, for all possible directed 
interactions (indicated by green or orange arrows). Boxes indicate the interquartile range and 
whiskers the confidence interval for the estimated TEs. TEs above the grey horizontal band 
indicate statistically significant causal influences (see Methods). In each plot, the third and the 
fourth boxes (from left to right) refer to TEs evaluated from "LFPs" restricted to groups of neu- 
rons that are source and target of a TL (pale green color denotes TL in the "green-to-orange" 
area direction, lilac color denotes TL in the "orange-to-green" area direction). Below each TE 
box-plot, effective connectivity is also represented in a diagrammatic form. Arrow thicknesses 
encode the strength of corresponding causal interactions (if statistically significant). D: TEs for 
the unidirectional driving effective motif with embedded TLs {Ktl = 22). E: TEs for the leaky 
driving effective motif with embedded TLs {Ktl = 24.5). Comparing these effective motifs 



37 



with Figures [3] and |4j we conclude that the embedding of TLs does not alter the overall effective 
connectivity. 

Supporting Figure S4. Scaling of Mutual Information (MI) with spike train length. 

MI normalized by entropy (at optimal time lag) is plotted against the inverse data fraction 
q. For each data fraction q, several bivariate spike trains are extracted from the original long 
spike trains (3 min, q = 1) and the mean MI is further averaged over these reduced-length 
spike trains. Asymptotic values are extrapolated through a quadratic interpolation. Error bars 
correspond to standard error. A: unidirectional driving effective motif, MI along the TL in the 
leader-to-laggard direction (pale green color), extrapolated asymptotic value is MI/H = 0.683. 
B: unidirectional driving effective motif, MI along the TL in the laggard-to- leader direction (lilac 
color), extrapolated asymptotic value is MI/H = 0.0066. In both cases, the finite size of the used 
spike trains produces a positive but small bias in the estimation of MI. Compared to Figure |8p, 
for the leader-to-laggard direction the overestimation is of ~ 3% and for the laggard-to-leader 
direction is of ~ 6%. 

Supporting Text SI. Full description of model parameters and complete analytic 
expressions. This text contains the following sections: i) Model neurons; ii) Model synapses; 
iii) Parameters of the background noise; iv) Phase response of the rate model; v) Phase-locking 
in the rate model. 
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Supporting Text SI 



Model neurons 



We use the Wang-Buzski (WB) conductance-based model (ref. [86] in main text) to describe 
each single excitatory and inhibitory neuron. The WB model is described by a single compart- 
ment endowed with sodium and potassium currents. The membrane potential is given by: 

C—j- = —II — iNa — Ik + lext + Irec 

at 

where C is the capacitance of the neuron, = giiy — ^l) is the leakage current, I^xt is 
an external driving current and Irec is due to recurrent interactions with other neurons in 
the network (see later). Sodiiim and potassium currents are voltage-dependent and given by 
iNa = ONa'mlohiV — Vno) and Ik = gK'n^iV — Vk)- The activation of the sodium current is 
instantaneous: 

m (V)~ 

Sodium current inactivation and potassium current activation evolve according to: 

- = ^.{a,{V){l-x)-^,{V)x) 

where x = h,n and ax and l3x{V) are non-linear functions of the membrane potential given by: 

0.1(F + 35) 



amiV) 



V+35 

1 + e 10 



, s V+60 

0.03(1/ + 34) 
1 — e 10 

, < V+44 

I3^(V) = 0.375e-^ 



3 



_ V+28 

1 + e 10 



Other parameters are gNa = 35 mS/cm^, V^a = 55 mV, gx = Q ms/cm^, Vk = —90 mV, 
gL = 0.1 mS/cm^, C = 1 /xF/cm^ and (f) = 5. 



Model synapses 

The synaptic current induced in a postsynaptic neuron by a single presynaptic action potential 
is given by Ispikeit) = ~gxSspike{t)iV — Vx), where V is the potential in the postsynaptic neuron 
and Vx is the reversal potential of the synapse (for excitatory synapses Ve = mV, for inhibitory 
synapses Vj = —80 mV). The time-course of the postsynaptic conductance is described by: 



Sspike{t) OC (exp {-{t + d- t*)/Ti) - exp {-{t + d- t*)/T2)) 
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for t > t*, otherwise, where t* is the time of the presynaptic spike, d is the latency, ri the 
rise-time and T2 the decay-time. The total recurrent current Irec{t) is the sum of time-dependent 
contributions Ispike{t) from all the presynaptic spikes fired to time t. The normalization constant 
of Sgpikeit) is chosen such as the peak value of Sgpike is equal to 1. For all simulations in the paper, 
we take n = 1 ms, r2 = 3 ms and d = 0.5 ms. Thus, post-synaptic currents have a relatively 
fast decay, corresponding to AMPA-like excitatory and GABA^i-like inhibitory synapses. For 
simplicity, we take only two possible peak conductances, gj = 90 /iS/cm^ for inhibitory synapses 
within an area and gE = ^ jjS/cvo? for excitatory synapses within and between areas. 

Parameters of the background noise 

In addition to recurrent synaptic inputs, each neuron receives a noisy input, representing back- 
ground spiking activity. It is modeled as an excitatory current having the same functional form 
of a recurrent current induced by a Poisson spike train with firing rate fext- The peak con- 
ductance of this noisy background input is gext- In our simulations, wc take f^xt = 5 kHz, 
and gext = 9e = 5 /iS/cm^. Each neuron is driven by statistically independent Poisson noise 
realizations. 

Phase response of the rate model 

As previously throughly reported in the Supplementary Material of ref. [42] (in main text), the 
firing rate of a single oscillating area (only local inhibitory coupling Kf < with delay D) can be 
derived analytically assuming that: (i) the total input current Itot{t) = I + KiR{t — D) is below 
threshold (i.e. negative) for a duration Tgt > D; (ii) the delay D and the oscillation period T 
fulfill the inequalities D < T — Tst < 2D. The conditions (i) and (ii) hold for sufficiently strong 
local inhibition, and, specifically, for the value Kj = —250 and the delay D = 0.1 adopted in the 
main paper. Under these conditions, the limit cycle of the firing rate assumes then the following 
analytic form (see Figure 2B in the main paper) : 
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where Rpeak is the peak ampHtude of the periodic oscillation of the rate and depends Hnearly 
on the level of the background current /. The oscillation period T and the sub-threshold time 
can be determined numerically by solving the system of non-linear equations: 
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We define the phase relative to the oscillation as 4>{t) = mod (t — to,T), where the time-shift to 
is chosen such as (t>{ttpeak) = in correspondence of the timings tpeak of oscillation peaks. Phases 
are therefore, with this notation, bounded between and 1. We use this convention throughout 
all analytic developments for the sake of simplicity. In The results involving phases in the main 
article are then translated back into the more usual angular range comprised between 0° and 
360°. The application of a pulse current 51 = h5{(t) — (pp) a, phase (^p induces a phase-shift 
S^{<pp) = hZ{(f)p) (see Figure S3B). The analytic expression for the Phase Response Curve (PRC) 
Z{(f)) can be derived from the knowledge of the limit cycle solution, and reads: 
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^ The resulting PRC is therefore null over a very large interval of 



where <j)st = -if- and ((>£) — ^ 
phases, leading in this broad range to refractoriness toward perturbations. A plot of Z 
the parameters used in our study is reported in Figure 4D (main text). 



for 



Phase-locking in the rate model 

As discussed in the main text, the time-evolution of the instantaneous phase shift A(f){t) between 
two coupled areas can be described, in the weak coupling limit, by the equation: 

The term r(A0) is a functional of the phase response and of the limit cycle waveform of the 
uncoupled oscillating areas. In terms of the previously derived analytic expressions of 2^(0) 
and of the rate oscillation limit cycle R{<J)) (phase-reduced) for Ke = 0, this functional can be 
expressed as r(A0) = C(A0) — C(— A0), where: 



Jo 



C{A(j))= / Z{(j))R{(f) + A(t)-D)d4> 
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Stable phase- lockings are therefore given by the zeroes of T with negative slope crossing. Analytic 
expressions for the integral C(A0) have already been derived and published in the Supplemen- 
tary Material of ref. [31] (in main text). We report here these expression again, in order to 
make the presentation of results self-contained. To compute C{A4>), six different intervals of 
Acf) need to be considered separately. The result is: 



C{Acj>) 



Cooifpst, 1) + Cloicpst, 1 - 4>d) 

Coo{(t>st, 1) + Cio{(t>st, l-(f>D) + Coii(t)st + (t>D- A0, 1) 



A(f) G [<^st + (f)D-l, <t>st + 2013 



Coo{<t)su 1) + Cio{<l)st, l-(t)D) + Coi{(^st + <Pd- Act), 1) 

+Cii{<j)st + (pD-A(f),l- <j)D) + Co2{'f>st + '^(l>D - A(f), 1) G [(f>st + - 1, M 

Coo{<Psu 1 + - A(/.) + e^Coo(0st + 0d - Ac/., 1) 
+Cio(0.t, 1-(^d) + Coi{(Pst, l + (t>D- A(P) 

+Cii(0,t, 1 - 0d) + Co2{(t)st + '^(t>D - A0, 1 + - A,^) Acj) e [<^d, + 



Cooicpst, l + (t>D-A(f>) + (^Cm{(pst + (^D- Act), 1) 
+Cio{ct)su 1-(^d) + CQi{ct)su l + (l>D- Act)) 
+Cii(0st, l-ct)D) + Co2{ct)st + 2<^D - Act),l + ct>D- Acf)) 

+Cl2(0st + 2(/.D - A0, 1 - ,^13) 

Coo(<^.t, 1 + - Ac/.) + e^Coo(0st + (/-D - Ac/., 1) 

+Cio(0st, 1 + c/.^ - Ac/.) + e^Cio(l + </>D - Ac/., 1 - c/.o) 
+^01 (</.,*, 1 + (/.D - Ac/.) + (:7ii(c/.,t, 1 + (/.D - Ac/.) 
+Co2(0st, 1 + (/"D - Ac/.) + Ci2(c/.st, 1 + </>!? - Ac/.) 



Ac/>G [c/.,t- l + 3c/>i3,2c/>i3] 



Ac/) > 
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where 

Coo (a, 6) 
Cio{a,b) 

Coi{a,b) 
Co2(a, 6) 



Cuia,b) 



Ci2ia,b) 



a)Te 



-T(l-A4>+<Pn) 



T{x + . 



If 



-Kje 



D-T 



_^gT(</,c-A</,) 



r(x + Ac/. - 2ct>D - ^st? ^ T{x + - 2<Pd - ^stf^ 



6 



T(0,-A<^) _ + - T){b - a)) 

_^gT(</,c-A</.) 



(13) 



^ + T{x + Ac/. - 2,/.^ - </..t)^ 



(14) 



(xr + g-r)3 (xr + zj-r)2 (xT + £>-r)^ 
+ r(i + A0 - - 05t) + h 



+r(l + Ac/. - 3<Pd - 0.t) ^^^^f + r(l + Ac/. - - 0.^)^ ^""^^? 



where [/(a;)]^ = /(6) — /(a). A plot of r(A0) for the parameters used in our study is reported 
in Figure 4B (main text). 
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Figure 3 
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Figure 4 
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Figure 8 




Figure 10. [FIGURE SI] 
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Figure 11. [FIGURE S2] 
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Figure 12. [FIGURE S3] 



1.0 
0.9 
0.8 
0.7 
0.6 



B 




0.02 r 

0.015 
0.01 

0.005 
0.00 



S5 



02468 024 
Inverse data fraction q 



Figure 13. [FIGURE S4] 



